library(vegan)
library(tidyverse)
library(dynutils)
otu <-
    read.table(
        file = "./metaphlan.tsv",
        sep = "\t",
        header = F,
        row.names = 1,
        quote = ""
    ) %>%
    t() %>%
    as_tibble()

group <-
    read.table(
        file = "./group_info.tsv",
        sep = "\t",
        header = T,
        # row.names = 1,
        quote = ""
    ) %>% filter(., !grepl("control", Group)) # remove control
# 导入临床数据(除了性别年龄)
clinic_original <-
    readxl::read_xlsx("./clincal.xlsx",
        sheet = 1
    ) %>%
    mutate(subject = str_sub(sample, 2, 3) %>% as.numeric() %>% as.character()) # 生成样本序号, 并转换为character类型
# import sex and age
clinic_sexage <-
    readxl::read_xlsx("./clincal.xlsx",
        sheet = 2
    ) %>%
    filter(!grepl("C.*", m0)) %>%
    pivot_longer(
        cols = 1:3,
        names_to = "Group",
        values_to = "samples",
        values_drop_na = T
    ) %>%
    select(-Group) %>%
    mutate(sex_factor = ifelse(sex == "male", 1, 2) %>% as.character()) # 将性别转换为character
clinic_original <-inner_join(clinic_original,clinic_sexage %>%  select(-sex),by=c("sample"="samples")) #合并临床指标和性别年龄
tax_s <-
    select(otu, matches(".*s__.*|clade")) %>%
    filter(clade_name %in% group$Sample) %>%
    rename(sample = clade_name)
# select genus from otu
tax_g <-
    select(otu, matches(".*g__((?!s__).)*$|clade", perl = T)) %>%
    filter(clade_name %in% group$Sample) %>%
    rename(sample = clade_name)
# select genus from otu
path <-
    read.table(
        file = "./path.csv",
        sep = ",",
        header = F,
        row.names = 1,
        quote = ""
    ) %>%
    t() %>%
    as_tibble() %>%
    filter(Pathway %in% group$Sample) %>%
    rename(sample = Pathway)
# --------------------------------------------------#
# 单个指标分别计算
# --------------------------------------------------#
single_xy <- function(dist.abund, y) {
    dist.temp <- dist(y %>% as.numeric(), method = "euclidean")
    abund <- mantel(dist.abund, dist.temp, method = "spearman", permutations = 9999, na.rm = TRUE)
    abund
}
get_single_mantel <- function(x, y) {
    # browser()
    x <- semi_join(x, y, by = ("sample"))
    y <- semi_join(y, x, by = ("sample"))
    dist.abund <- vegdist(x %>% select(-c("sample")) %>% transmute_all(as.numeric), method = "bray")
    y_list <- as.list(y %>% select(-c("sample")))
    result <- map2(rep(list(dist.abund), length(y_list)), y_list, single_xy) %>% set_names(names(y_list))
    result
}
single_result <- map2(list(tax_s, tax_g, path), list(clinic_original, clinic_original, clinic_original), get_single_mantel) %>% set_names(c("tax_s", "tax_g", "path"))
# 多种变量协同作用
get_multi_mantel <- function(x, y) {
    # browser()
    x <- semi_join(x, y, by = ("sample"))
    y <- semi_join(y, x, by = ("sample"))
    dist.abund <- vegdist(x %>% select(-sample) %>% transmute_all(as.numeric), method = "bray")
    dist.temp <- dist(y %>% select(-sample, -subject,-sex_factor) %>% scale(center = TRUE, scale = TRUE)
        %>% as_tibble(), method = "euclidean")
    abund <- mantel(dist.abund, dist.temp, method = "spearman", permutations = 9999, na.rm = TRUE)
    abund
}
multi_result <- map2(list(tax_s, tax_g, path), list(clinic_original, clinic_original, clinic_original), get_multi_mantel) %>% set_names(c("tax_s", "tax_g", "path"))
# 整理结果到数据框
to_tibble <- function(single) {
    tmp <- map(single, ~ tibble(r = .$statistic, p = .$signif)) %>% set_names(names(single))
    tib <- tibble(r = map(tmp, "r") %>% unlist(), p = map(tmp, "p") %>% unlist()) #  %>%  mutate(R=map(r,))
    tib
}
all_result <- map(single_result, to_tibble)
# 合并结果
all_result[[1]] <- add_row(all_result[[1]],
        r=c(multi_result[[1]][["statistic"]]),
        p=c(multi_result[[1]][["signif"]]))
all_result[[2]] <- add_row(all_result[[2]],
        r=c(multi_result[[2]][["statistic"]]),
        p=c(multi_result[[2]][["signif"]]))
all_result[[3]] <- add_row(all_result[[3]],
        r=c(multi_result[[3]][["statistic"]]),
        p=c(multi_result[[3]][["signif"]]))
write.table(all_result,"./mantel.csv", sep=",", quote=F)
